Figures for Seasonally variable thermal performance curves prevent adverse effects of heatwaves

Main Text Figures

comb_epr_plot = comb_preds %>% 
  filter(metric == "EPR") %>% 
  ggplot() +
  geom_point(aes(temp, rate, colour = curve_id), 
             filter(comb_d, metric == "EPR"), 
             size = 1.5, alpha = 0.6, 
             position = position_jitter(width = 0.5, height = 0)) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper, group = curve_id), 
              filter(comb_boot_conf_preds, metric == "EPR"), 
              fill = 'grey60', alpha = 0.3) +
  geom_line(aes(temp, .fitted, col = curve_id), linewidth = 2) +
  scale_colour_manual(values = comb_colors) + 
  labs(x = "", 
       y = "Egg Production \n(eggs/female/day)",
       colour = "Month") + 
  theme_matt(base_size = 12)

comb_hs_plot = comb_preds %>% 
  filter(metric == "HF") %>% 
  ggplot() +
  geom_point(aes(temp, rate, colour = curve_id), filter(comb_d, metric == "HF"), size = 1.5, alpha = 0.6, 
             position = position_jitter(width = 0.5, height = 0)) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper, group = curve_id), filter(comb_boot_conf_preds, metric == "HF"), fill = 'grey60', alpha = 0.3) +
  geom_line(aes(temp, .fitted, col = curve_id), linewidth = 2) +
  #scale_color_brewer(type = "div", palette = 5, direction = -1) + 
  #scale_color_viridis_d(option = "mako") + 
  scale_colour_manual(values = comb_colors) + 
  labs(x = "", 
       y = "Hatching Success \n(%)",
       colour = "Month") + 
  theme_matt(base_size = 12)

comb_rf_plot = comb_preds %>% 
  filter(metric == "RF") %>% 
  ggplot() +
  geom_point(aes(temp, rate, colour = curve_id), filter(comb_d, metric == "RF"), size = 1.5, alpha = 0.6, 
             position = position_jitter(width = 0.5, height = 0)) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper, group = curve_id), filter(comb_boot_conf_preds, metric == "RF"), fill = 'grey60', alpha = 0.3) +
  geom_line(aes(temp, .fitted, col = curve_id), linewidth = 2) +
  #scale_color_brewer(type = "div", palette = 5, direction = -1) + 
  #scale_color_viridis_d(option = "mako") + 
  scale_colour_manual(values = comb_colors) + 
  labs(x = "Temperature (°C)", 
       y = "Offspring Production \n(offspring/female/day)",
       colour = "Month") + 
  theme_matt(base_size = 12)

comb_tsc = ggplot(comb_surv, aes(x=Temp, y=Surv, colour=Month)) + 
  geom_point(size=1.5, position=position_jitter(width=0.1, height=0.03)) +
  xlab("Temperature (°C)")+
  ylab("Survivorship \n(proportion survived)")+
  labs(colour = "Month") + 
  geom_hline(yintercept = 0.5, linetype = "dashed") +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=T, linewidth = 2) +
  scale_y_continuous(breaks = c(0,1)) + 
  #scale_color_brewer(type = "div", palette = 5, direction = -1) + 
  #scale_color_viridis_d(option = "mako") + 
  scale_colour_manual(values = comb_colors) + 
  theme_matt(base_size = 12)

ggarrange(comb_epr_plot, comb_hs_plot, comb_rf_plot, comb_tsc, 
          ncol = 2, nrow = 2,
          common.legend = T, legend = "bottom", 
          labels = "AUTO", vjust = 1)

# ggarrange(comb_epr_plot, comb_hs_plot, comb_rf_plot, comb_tsc, nrow = 1,
#           common.legend = T, legend = "bottom")

combined_opt_coll = comb_params %>% 
  filter(metric == "RF" & term == "topt") %>% 
  ggplot(aes(x = growth_temp, y = estimate, shape = species)) + 
  geom_smooth(data = filter(comb_params, metric == "RF" & term == "topt" & curve_id != "November_2"),
              method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  scale_shape_manual(values = c(16,21)) + 
  ylab("Optimum (°C)") + 
  xlab("Collection Temperature (°C)") + 
  labs(colour = "Month") + 
  theme_matt(base_size = 12)

combined_opt_diff = comb_params %>% 
  filter(metric == "RF" & term == "topt") %>% 
  ggplot(aes(x = growth_temp, y = margin, shape = species)) +
  geom_hline(yintercept = 0, linewidth =1, linetype = "dashed") +
  geom_smooth(data = filter(comb_params, metric == "RF" & term == "topt" & curve_id != "November_2"),
              method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  ylab("Margin (°C)") + 
  xlab("Collection Temperature (°C)") + 
  scale_shape_manual(values = c(16,21)) + 
  theme_matt(base_size = 12) 

combined_ld_coll = ggplot(combined_tolerance, aes(x = Coll_temp, y = LD50, shape = species)) + 
  geom_smooth(method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  scale_shape_manual(values = c(16,21)) + 
  ylab("Thermal Tolerance (°C)") + 
  xlab("Collection Temperature (°C)") + 
  theme_matt(base_size = 12)

combined_ld_diff = ggplot(combined_tolerance, aes(x = Coll_temp, y = margin, shape = species)) +
  geom_smooth(method = "lm", colour = "grey50") + 
  geom_point(size = 3, stroke = 1) + 
  scale_shape_manual(values = c(16,21)) + 
  ylab("Warming Tolerance (°C)") + 
  xlab("Collection Temperature (°C)") + 
  theme_matt(base_size = 12)

ggarrange(combined_opt_coll, combined_opt_diff, combined_ld_coll, combined_ld_diff, ncol = 2, nrow = 2, common.legend = T,
          legend = "bottom", labels = "AUTO")

x.axis_labels = c("1" = "short", "2" = "long", "3" = "short", "4" = "long")

F0_grid = F0_rf_summary %>%
  mutate(month = fct_relevel(month, "June", "August", "November")) %>%
  ggplot(aes(x = duration, y = difference, colour = trait, shape = duration)) +
  facet_grid(. ~ month) +
  geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
  geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1.3) +
  geom_point(size = 5, fill = "white") +
  scale_colour_manual(values = c("body size" = "grey75", "production" = "black")) +
  scale_shape_manual(values = c("long" = 16, "short" = 21)) +
  scale_x_discrete(labels= x.axis_labels) +
  ggtitle("Direct Effects (F0)") +
  xlab("") +
  ylab("Effect Size\nHeatwave - Control") +
  ylim(-1,1.1) +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5),
        legend.position = "none")

F1_summary = bind_rows(F1_rf_effect_size, F1_bs_effect_size) %>%
  dplyr::select(variable, difference,
                bca_ci_low, bca_ci_high,
                month, duration, trait, generation, off_temp) %>%
  mutate("order_code" = paste(trait, duration, sep = "_"),
         "order_number" = case_when(
           order_code == "production_short" ~ 1,
           order_code == "production_long" ~ 2,
           order_code == "body size_short" ~ 3,
           order_code == "body size_long" ~ 4),
         month = fct_relevel(month, "June", "August", "November"))

F1_summary$order_number = factor(F1_summary$order_number, levels = c(1,2,3,4))
F1_grid = ggplot(F1_summary, aes(x = order_number, y = difference, colour = trait, shape = duration, group = trait)) +
  facet_grid(off_temp ~ month, ) +
  geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
  geom_line() +
  geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1) +
  geom_point(size = 3, fill = "white") +
  scale_colour_manual(values = c("body size" = "grey75", "production" = "black")) +
  scale_shape_manual(values = c("long" = 16, "short" = 21)) +
  xlim(0.5,4.5) +
  scale_x_discrete(labels= x.axis_labels) +
  xlab("") +
  ylab("Effect Size\nHeatwave - Control") +
  ggtitle("Transgenerational Effects (F1)") +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank(),
        strip.background.x = element_blank(),
        strip.text.x = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5))

ggarrange(F0_grid, F1_grid, nrow = 2, ncol = 1, heights = c(0.45,1), common.legend = T, legend = "right")

Supplemental Information

#field tpc parameters
comb_params %>%  
  mutate(curve_id = fct_relevel(curve_id, c("January", "February", "March", "April", "May", "June", 
                                            "July", "August", "September", "October", "November_1", "November_2"))) %>% 
  ggplot(aes(x = curve_id, y = estimate, colour = species)) +
  facet_wrap(term~metric, scales = 'free_y') + 
  geom_point(size = 4) +
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  scale_colour_manual(values = c("royalblue1", "indianred2")) + 
  labs(x = "Month",
       y = "Parameter Estimate",
       colour = "Species") + 
  theme_bw(base_size = 16) +
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 315, hjust = 0, vjust = 0.5))

#field tpc parameters
ggplot(comb_params, aes(x = growth_temp, y = estimate, colour = species)) +
  facet_wrap(term~metric, scales = 'free_y') + 
  geom_smooth(data = filter(comb_params, curve_id != "November_2"),
              method = "lm", se = F) + 
  geom_point(size = 4) +
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  scale_colour_manual(values = c("royalblue1", "indianred2")) + 
  labs(x = "Collection Temperature (°C)",
       y = "Parameter Estimate",
       colour = "Species") + 
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank())

f1_size_data = F1_fbs %>% 
  ungroup() %>% 
  mutate(Day = if_else(Day == "1_to_3", "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November"),
         Day = fct_relevel(Day, "Short", "Long"),
         Parental_treatment = fct_relevel(Parental_treatment, "Heatwave", "Control"))


F1_size.model = lm(data = f1_size_data,
                    Size ~ Parental_treatment * Day * Offspring_temp * Month)

performance::check_model(F1_size.model, check = c("pp_check", "linearity", "homogeneity", "outliers", "qq"))

tsm = comb_params %>% 
  filter(metric == "RF" & term == "topt") %>% 
  filter(margin < 15) %>% 
  select("Month" = curve_id, species, "Coll_temp" = growth_temp, margin)

wt = combined_tolerance %>% 
  select(Month, species, Coll_temp, margin)

tsm_slopes = lm(data = tsm, margin ~ Coll_temp * species) %>% 
  emmeans::emtrends(~species, var = "Coll_temp") %>% 
  data.frame()

tsm_slope_plot = ggplot(tsm_slopes, aes(x = species, y = Coll_temp.trend)) + 
  geom_hline(yintercept = 0) + 
  geom_point(size = 4) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),  
                width = 0.1, linewidth = 1) + 
  labs(x = "Species",
       y = "Safety Margin Slope\n(°C/°C)") + 
  theme_matt()

wt_slopes = lm(data = wt, margin ~ Coll_temp * species) %>% 
  emmeans::emtrends(~species, var = "Coll_temp") %>% 
  data.frame()

wt_slope_plot = ggplot(wt_slopes, aes(x = species, y = Coll_temp.trend)) + 
  geom_hline(yintercept = 0) + 
  geom_point(size = 4) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),  
                width = 0.1, linewidth = 1) + 
  labs(x = "Species",
       y = "Warming Tolerance Slope\n(°C/°C)") + 
  theme_matt()

ggarrange(tsm_slope_plot, wt_slope_plot, labels = "AUTO")

f0_model_females = F0_epr %>% 
  group_by(Month, Treatment, Female) %>% 
  count() %>% 
  filter(n == 2) %>% 
  mutate('female_ID' = paste(Month, Treatment, Female, sep = "_"))

f0_model_data = F0_epr %>% 
  mutate('female_ID' = paste(Month, Treatment, Female, sep = "_")) %>% 
  filter(female_ID %in% f0_model_females$female_ID) %>% 
  mutate(Day = if_else(Day == "1_to_3", "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November"),
         Day = fct_relevel(Day, "Short", "Long"))

f1_model_data = F1_epr %>% 
  mutate(Offspring_temp = as.factor(Offspring_temp)) %>% 
  ungroup() %>% 
  mutate(Day = if_else(Day == "1_to_3", "Short", "Long"),
         Month = fct_relevel(Month, "June", "August", "November"),
         Day = fct_relevel(Day, "Short", "Long")) %>% 
  mutate(Parental_treatment = fct_relevel(Parental_treatment, "Heatwave", "Control"))

f0_epr_plot = ggplot(f0_model_data, 
                 aes(x = Day, y = Total, colour = Treatment)) + 
  facet_grid(.~Month, scales = "free_y") + 
  geom_hline(yintercept = 0) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  labs(y = "Egg Production \n(per female per day)") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        axis.title.x = element_blank(),
        strip.text = element_text(size = 18))

f1_epr_plot = ggplot(f1_model_data, 
                 aes(x = factor(Offspring_temp), y = Total, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  #geom_violin(position = position_dodge(width = 1)) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  scale_y_continuous(breaks = c(0, 100, 200)) + 
  labs(x = "Offspring Temp. (°C)", 
       y = "Egg Production \n(per female per day)",
       colour = "Parental Treatment") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        strip.text = element_text(size = 18))

raw_epr = ggarrange(f0_epr_plot, f1_epr_plot, nrow = 2, 
          heights = c(0.3, 0.7), common.legend = T, legend = "bottom")

###

f0_hs_plot = ggplot(f0_model_data, 
                 aes(x = Day, y = Success, colour = Treatment)) + 
  facet_grid(.~Month, scales = "free_y") + 
  geom_hline(yintercept = 0) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  labs(y = "Hatching Success") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        axis.title.x = element_blank(),
        strip.text = element_text(size = 18))

f1_hs_plot = ggplot(f1_model_data, 
                 aes(x = factor(Offspring_temp), y = Success, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  #geom_violin(position = position_dodge(width = 1)) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  scale_y_continuous(breaks = c(0, 1)) + 
  labs(x = "Offspring Temp. (°C)", 
       y = "Hatching Success",
       colour = "Parental Treatment") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        strip.text = element_text(size = 18))

raw_hs = ggarrange(f0_hs_plot, f1_hs_plot, nrow = 2, 
          heights = c(0.3, 0.7), common.legend = T, legend = "bottom")

###

f0_prod_plot = ggplot(f0_model_data, 
                 aes(x = Day, y = Hatched, colour = Treatment)) + 
  facet_grid(.~Month, scales = "free_y") + 
  geom_hline(yintercept = 0) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  labs(y = "Offspring Production \n(per female per day)") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        axis.title.x = element_blank(),
        strip.text = element_text(size = 18))

f1_prod_plot = ggplot(f1_model_data, 
                 aes(x = factor(Offspring_temp), y = Hatched, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  #geom_violin(position = position_dodge(width = 1)) + 
  geom_boxplot(position = position_dodge(width = 0.5),
               width = 0.3) + 
  geom_point(position = position_dodge(width = 0.5),
             alpha = 0.5) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  scale_y_continuous(breaks = c(0, 100, 200)) + 
  labs(x = "Offspring Temp. (°C)", 
       y = "Offspring Production \n(per female per day)",
       colour = "Parental Treatment") + 
  theme_matt() + 
  theme(panel.border = element_rect(fill = NA, colour = "black"),
        legend.position = "bottom",
        strip.text = element_text(size = 18))

raw_production = ggarrange(f0_prod_plot, f1_prod_plot, nrow = 2, 
          heights = c(0.3, 0.7), common.legend = T, legend = "bottom")

ggarrange(raw_epr, raw_hs, raw_production, nrow = 1, common.legend = T, legend = "bottom",
          labels = "AUTO")

#Subsequent Rows
F1_hs_effect_size$trait = "HS"
F1_hs_effect_size$generation = "F1"

F1_total_effect_size$trait = "EPR"
F1_total_effect_size$generation = "F1"

F1_rf_effect_size$trait = "OP"
F1_rf_effect_size$generation = "F1"

F1_bs_effect_size$trait = "Size"
F1_bs_effect_size$generation = "F1"

F0_size_summary = data.frame("trait" = "Size", generation = "F0", month = "June", duration = "long")

F0_data = bind_rows(F0_hs_summary, F0_total_summary,F0_rf_summary, F0_size_summary) %>%
  dplyr::select(trait, difference, bca_ci_low, bca_ci_high, month, duration, trait, generation) %>%
  mutate("order_code" = paste(trait, duration, sep = "_"),
         "order_number" = case_when(
           order_code == "total_short" ~ 1,
           order_code == "total_long" ~ 2,
           order_code == "success_short" ~ 3,
           order_code == "success_long" ~ 4,
           order_code == "production_short" ~ 5,
           order_code == "production_long" ~ 6,
           order_code == "size_long" ~ 7),
         month = fct_relevel(month, "June", "August", "November"),
         trait = case_when(
           trait == "total" ~ "EPR", 
           trait == "success" ~ "HS", 
           trait == "production" ~ "OP", 
           trait == "Size" ~ "Size"),
           fct_relevel(trait, "EPR", "HS", "OP", "Size"),
         duration = fct_relevel(duration, "short", "long"),
         group_ID = paste(month, trait, sep = "_"))

F0_data$order_number = factor(F0_data$order_number, levels = c(1,2,3,4,5,6,7))


F1_data = bind_rows(F1_total_effect_size, F1_hs_effect_size, F1_rf_effect_size, F1_bs_effect_size) %>%
  dplyr::select(trait, difference, bca_ci_low, bca_ci_high, month, duration, generation, off_temp) %>%
  mutate("order_code" = paste(trait, duration, sep = "_"),
         "order_number" = case_when(
           order_code == "EPR_short" ~ 1,
           order_code == "EPR_long" ~ 2,
           order_code == "HS_short" ~ 3,
           order_code == "HS_long" ~ 4,
           order_code == "OP_short" ~ 5,
           order_code == "OP_long" ~ 6,
           order_code == "Size_short" ~ 7,
           order_code == "Size_long" ~ 8),
         trait = if_else(trait == "epr", "total", trait),
         month = fct_relevel(month, "June", "August", "November"),
         trait = fct_relevel(trait, "EPR", "HS", "OP", "Size"),
         duration = fct_relevel(duration, "short", "long"))


F1_data$order_number = factor(F1_data$order_number, levels = c(1,2,3,4,5,6,7,8))

#Top row - F0 (direct effects)
x.axis_labels = c("1" = "short", "2" = "long", "3" = "short", "4" = "long",
                  "5" = "short", "6" = "long", "7" = "short", "8" = "long")

F0_grid = ggplot(F0_data, aes(x = duration, y = difference, colour = trait, group = group_ID)) +
  facet_grid(. ~ month) +
  geom_line(position = position_dodge(width = 0.7),
            linewidth = 1) +
  geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
  geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1,
                position = position_dodge(width = 0.7)) +
  geom_point(size = 4, fill = "white", position = position_dodge(width = 0.7)) +
  scale_shape_manual(values = c("long" = 16, "short" = 21)) +
  scale_x_discrete(labels= x.axis_labels) +
  scale_colour_manual(values = c("Size" = "darkgrey",
                                 "HS" = "gold",
                                 "OP" = "forestgreen",
                                 "EPR" = "cornflowerblue")) +
  xlab("") +
  ylab("Effect Size\nHeatwave - Control") +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5),
        legend.position = "none")

#Following three rows - F1 (transgeneration / indirect effects)
F1_grid = ggplot(F1_data, aes(x = duration, y = difference, colour = trait, group = trait)) +
  facet_grid(off_temp ~ month, ) +
  geom_hline(yintercept = 0, colour = "black", linewidth = 0.3) +
  geom_line(position = position_dodge(width = 0.5),
            linewidth = 1) +
  geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high), width = 0, linewidth = 1,
                position = position_dodge(width = 0.5)) +
  geom_point(size = 3, fill = "white", position = position_dodge(width = 0.5)) +
  scale_shape_manual(values = c("long" = 16, "short" = 21)) +
  xlim(0.5,4.5) +
  scale_x_discrete(labels= x.axis_labels) +
  scale_colour_manual(values = c("Size" = "darkgrey",
                                 "HS" = "gold",
                                 "OP" = "forestgreen",
                                 "EPR" = "cornflowerblue")) +
  xlab("") +
  ylab("Effect Size \nHeatwave - Control") +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank(),
        strip.background.x = element_blank(),
        strip.text.x = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(angle = -45, hjust = 0, vjust = 0.5))

ggarrange(F0_grid, F1_grid, nrow = 2, ncol = 1, heights = c(0.35,1), common.legend = T, legend = "right", labels = "AUTO")

F0_rf_summary$month = factor(F0_rf_summary$month, levels = c("June", "August", "November"))
F0_rf_summary$duration = factor(F0_rf_summary$duration, levels = c("short", "long"))
F0_dur_summary$month = factor(F0_dur_summary$month, levels = c("June", "August", "November"))

ggplot(F0_dur_summary, aes(x = month, y = difference, fill = treatment)) + 
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity", 
           position = position_dodge(width = 0.9),
           colour = "black") + 
  geom_errorbar(aes(ymin = bca_ci_low, ymax = bca_ci_high),
                position = position_dodge(width = 0.9),
                width = 0.2, linewidth = 0.75) + 
  scale_fill_manual(values = c("Control" = "grey30", "Heatwave" = "white")) + 
  labs(y = "Effect Size (Hedge's g)\nLong - Short events",
       x = "") + 
  theme_matt()

effect_corr = F1_summary %>%
  select(trait, difference, month, duration, off_temp) %>%
  pivot_wider(id_cols = c(month, duration, off_temp),
              names_from = trait,
              values_from = difference)

ggplot(effect_corr, aes(x = `body size`, y = production)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = F,
              colour = "grey60",
              size = 1) +
  labs(x = "Body Size Effect",
       y = "Production Effect") +
  theme_matt()

day_cols = c("Short" = "grey80", "Long" = "grey30")

size_temp1 = ggplot(f1_size_data, aes(x = Offspring_temp, y = Size, colour = Parental_treatment)) + 
  facet_grid(Month~Day) + 
  geom_jitter(width = 0.5, size = 1.6, alpha = 0.4) +
  geom_smooth(method = "lm", linewidth = 1.4, alpha = 0.2) + 
  labs(x = "Offspring Temperature (°C)",
       y = "Size (mm)",
       colour = "Parental Treatment") + 
  scale_x_continuous(breaks = c(12,17,22)) + 
  scale_colour_manual(values = c("Heatwave" = "coral3", "Control" = "skyblue3")) + 
  theme_matt(base_size = 15) + theme(legend.position = "bottom",
                                     panel.grid = element_blank(),
                                     panel.border = element_rect(fill = NA, colour = "black"))

size_temp2 = emmeans::emtrends(F1_size.model, c("Month", "Day", "Parental_treatment"), var = "Offspring_temp") %>% 
  as.data.frame() %>% 
  ggplot(aes(x = Parental_treatment, y = Offspring_temp.trend, 
             colour = Day, group = Day)) + 
  facet_wrap(Month~.) + 
  geom_hline(yintercept = 0) + 
  geom_line(linewidth = 1.5,
            position = position_dodge(width = 0.5)) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                linewidth = 1.5, width = 0.5,
                position = position_dodge(width = 0.5)) + 
  geom_point(size = 3,
             position = position_dodge(width = 0.5)) + 
  scale_colour_manual(values = day_cols) + 
  labs(x = "Parental Treatment",  
       y = "Size Slope (mm / °C)") + 
  theme_matt(base_size = 15) + 
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 290,
                                   hjust = 0, vjust = 0.5),
        plot.background = element_rect(fill = "white"))

ggarrange(size_temp1, size_temp2, labels = "AUTO", common.legend = F, legend = "bottom")

f1_size_contrasts = emmeans::emmeans(F1_size.model, ~ Parental_treatment | Day * Month * Offspring_temp, type = "response", at = list(Offspring_temp = c(12,17,22))) %>% 
  pairs() %>% data.frame() %>% 
  drop_na() %>% 
  mutate("ID" = paste0(Month, Offspring_temp)) %>% 
  filter(!(ID == "November12" & Day == "Short")) %>% 
  filter(!(ID == "November22" & Day == "Short"))

ggplot(f1_size_contrasts, aes(x = Day, y = estimate, group = ID)) +
  facet_grid(Offspring_temp~Month) + 
  geom_hline(yintercept = 0) + 
  geom_line(linewidth = 0.8) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE), width = 0.1, linewidth = 1) + 
  labs(x = "Duration",
       y = "Effect (mm; Control - Heatwave)") + 
  theme_bw() + 
  theme(panel.grid = element_blank())